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Abstract 

We present molecular-dynamic simulations of memory resistors (memristors) including the crys- 
tal field effects on mobile ionic species such as oxygen vacancies appearing during operation of the 
device. Vacancy distributions show different patterns depending on the ratio of a spatial period of 
the crystal field to a characteristic radius of the vacancy- vacancy interaction. There are signatures 
of the orientational order and of spatial voids in the vacancy distributions for some crystal field 
potentials. The crystal field stabilizes the patterns after they are formed, resulting in a non-volatile 
switching of the simulated devices. 

PACS numbers: 71.38.-k, 74.40.+k, 72.15.Jf, 74.72.-h, 74.25.Fy 
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Electrical switching in thin-film oxide nanodevices attracts renewed attention potentially 
enabling scaling of logic and memory circuits well beyond the current state of the art tech- 
niques [lj. The microscopic nature of resistance switching and a charge transport in such 
devices is still a subject of debate, but there is growing evidence that nonvolatile switching 
requires some sort of atomic rearrangement of mobile ionic species that drastically affects 
the current. The microscopic understanding of the atomic-scale mechanism and identifica- 
tion of the material changes within the device appears to be invaluable for controlling and 
improving the memristor performance. 

The number of oxygen vacancies within the volume 10x10x2 nm 3 of perspective nano- 
memristors, such as based on an amorphous layer of titanium dioxide, TiC^-x, could be 
as small as about a thousand, so that the conventional statistical (diffusion) approach for 
dealing with many-particle systems may fail. Recently, we have proposed a model for the 
kinetic behavior of oxide memristors and studied it using the Molecular Dynamics (MD) 
simulations of the Langevin equations describing the thermal diffusion and drift of inter- 
acting oxygen vacancies 2|. Our MD simulations revealed a significant departure of the 
vacancy distributions across the device from that expected within a standard drift-diffusion 
approximation. The channel formation in systems like TiC>2, NiO, VO2 is certainly accom- 
panied by local heating that is witnessed by the local emergence of high-temperature phases 
and observed by thermal microscopy and other techniques . Accounting for local heat- 
in g , we extended our MD modeling of oxygen vacancies driven by an externa, bias voltage 
by taking into account temperature gradients in thin films of oxide memristors [2J. Tem- 
perature gradients, producing strong variations of the rate of diffusion, affect the vacancy 
patterns in memristors. Our simulations indicated that variations of temperature across 
the memristor favor the formation of short-circuiting or shunting vacancy channels. At the 
same time, considerable temperature gradients along the sample can, by contrast, produce 
vacancy-poor regions where vacancy shunts are not formed and electron percolation paths 
are blocked. Here we extend our MD simulations of the memristor including the crystal field 
in the Langevin equations. 

Observations of oxygen vacancy migration and clustering in bulk [6j and nanoscale 3|, [7] 
samples of Ti02, induced by an electric- field, allows us to model a memristor [2] with the 
vacancy interaction potential, Fig. (TJ with the results shown in Fig. [2j In the model, there 
is a reduced rutile thin layer Ti0 2 - x near one of the metallic electrodes stabilized by the 
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FIG. 1: (Color online) The Lennard-Jones potential (steeper curve) compared with the O 2 
vacancy- vacancy interaction potential, Eq.(|2]). 

Coulomb mirror potential. Vacancies from this layer can drift toward the opposite electrode 
pushed by a pulse of an electric field. Vacancies interact with each other via the pairwise 
potential W, and with the electric field corresponding to the time-dependent deterministic 
force F. The environment exerts an independent Gaussian random force, £ on each particle 
with zero mean and intensity controlled by the temperature. Different from our previous 
studies (2I, we now include the periodic crystal field potential Z7(aCj) in the overdamped 
Langevin equations describing the drift-diffusion of the i-th particle as 

V — - F i ( Xi , t) - ^ - ^ — ^ + y/2k B Tr& (t). (1) 

where x" is the a-coordinate of the zth vacancy (a = x,y,z), rj is the friction coefficient, F-* 
is the a-component of electric pulse force, and is a random force. 

The 2 ~ vacancy-vacancy interaction potential, which is responsible for the clustering, 
is often modeled as [8j 

W{R) = Aexp(-R/a ) - B/R 6 + e 2 /(7re e J R), (2) 



where the short-range repulsive and attractive parts are represented with the parameters 



A = 22764.0 eV, a = 0.01490 nm, and B = 27.88 x 10" 6 eV nm 6 in Ti0 2 _x and the long- 
range Coulomb repulsion is given by the last term, see FigJTJ Since we simulate a limited 
number of vacancies, one can refer to each particle in our simulations as a cluster of vacancies, 
where a cluster-cluster interaction is more conveniently described by the combination of the 
Lennard- Jones and Coulomb potentials acting with the force 

F(R) = {^E L j[(R min / R) 12 — (R min /R) 6 ] + E c R min /R] , (3) 

where the relative strength of the Coulomb potential is given by E c / Elj = 2. This results in 
the position of the potential maxima R ma x ~ ^Rmin and the height of the potential barrier 
on the order of the depth of the potential well. The electric pulse strength is taken such that 
the Coulomb force is about ten times stronger than the maximum attracting force between 
vacancies, and the interactions are cut-off at a distance of about R m i n /20. The Lennard- 
Jones potential provides somewhat stiffer repulsion at small distances, but otherwise it fairly 
fits Eq.(J2J) with the parameters E LJ « 2.16 eV and R m in ~ 0.125 nm, FigJTJ 

The vacancy distribution depends on the boundary conditions, the sample size, tempera- 
ture gradients j^J, and, in the present simulations, on the ratio of the lattice constant, a and 
the characteristic scale of the interaction potential Rmin- Hence, it is most convenient to 
measure coordinates in units of R m %n by introducing r = x/i? m j n while we measure time in 
units of t = R m inl D by introducing r = t/t , where D = k B T/i] is the diffusion coefficient. 
The Langevin equations in this dimensionless (r, r) space takes the following form: 

drf d[V(r t ) + U(r t )}/(k B T) dW( ri - rj )/(k B T) ^ 

= &a 2. frf + (r), (4) 

where V(ri) is the electric-field potential, and the components of the (dimensionless) random 
force, satisfies the fluctuation - dissipation relation (£f (0)£^ (r)) = S(T)5 a> p8ij. While the 
diffusion in Ti02 is anisotropic, for simplicity sake here and below we consider the uniform 
temperature across the device and the isotropic diffusion of vacancies in a two dimensional 
(2D) square lattice with the crystal field, 

U(r) = U sin(2vrx//) sin(2vn//Z), (5) 

where x, y are two components of r, and I = a/R min . 
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FIG. 2: (Color online) Vacancy dynamics in different crystal fields (upper panels (a): I = 0.2, lower 
panels (b): I = 2) at different times after the electric pulse has been applied: t = tq [a(l),b(l)], 
r = 2r [a(2),b(2)], r = 3r [a(3),b(3)], and r = 4r [a(4),b(4)]. 

We simulate Eqs.(|4]) for N = 70 vacancies with I = 2, U — 4E LJ /3 and with I = 0.2, 
Uq = E LJ /6, respectively, so that the forces (~ U /a) are of the same order of magnitude 
in both cases. We place all the vacancies randomly near the bottom of a toy sample and 
then let them evolve according to Eq. inside a rectangular box, which mimics the actual 
sample. We use the 2D simulation area with the aspect ratio L y /L x = 2 and periodic 
boundary conditions (BC) along the x— direction. Note that the periodic BC allow us to 
simulate an infinite area sample using a rather small number of particles. To use periodic 
BC, we include periodic images of the vacancies with respect to vertical boundaries of the 
simulation box. We also incorporate opposite polarity charges by adding mirror images of 
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vacancies with respect to the top and the bottom of the sample. 

The results of simulations with the homogeneous temperature [ksT « 1/30 of the in- 
teraction potential well, FigJT], are presented in Fig. [2J Different columns refer to different 
moments in time [time intervals are shown relative to the rectangular electric pulse dura- 
tion, tq]. The magnitude of the dimensionless electric pulse force is taken as L y /2r , so that 
vacancies are pushed by the pulse to the center of the sample. The distribution is spread 
out and its maximum moves towards the top of the sample with a constant velocity. 

The important effect of the crystal field is that clusters appear more ordered than those 
without the crystal field Q]. Namely, there are signatures of an orientational order at the 
angle of 45° for I = 0.2, and of spatial voids in the vacancy distribution on the order of 
a few lattice constants for 1 = 2. This indicates that a large variety of different vacancy 
phases (from liquid-like via smectic or nematic liquid crystal-like to crystal-like phases) can 
be achieved in proper domains of parameters (such as temperature, vacancy concentration, 
vacancy pinning, electric pulse intensity, sample shape etc.) in this strongly interacting 
system with competing thermal and quenched disorders. Actually, our MD simulations 
may shed a light on what kind of microscopic distributions could be behind a so-called 
"random circuit breaker network" behavior phenomenologically introduced in Ref. [h]] used 
to explain both bipolar and unipolar resistance switching in the Pt/SrTiO x /Pt capacitor 
and changeover between the two regimes. 

Here, it is worth mentioning an analogy between vacancies in memristors and supercon- 
ducting vortices, since in both cases long-range interactions, which can be both repulsive 
and attractive, and thermal/quenched disorder play a crucial role. It is well known 11] that 
a huge variety of vortex phases exists in superconductors. In addition to the standard vortex 
crystal and vortex liquid phases, different types of smectic phases were predicted and ob- 
served in superconductors. For anisotropic superconductors, where the interaction between 
vortices can change from repulsive to attractive at short distances, vortex chains were pre- 
dicted and observed that are reminiscent of the simulated filament-like clusters of vacancies 
in memristors. Moreover, vortex avalanches can form due to a competition between ther- 
mal/quenched disorder and vortex- vortex interactions that again has a deep analogy with 
the observed vacancy clusters in memristors. 

Finally, it is well known that in addition to the equilibrium vortex phases, many different 
metastable phases can occur when vortices are driven by a magnetic field and relaxation 
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FIG. 3: (Color online) Vacancy dynamics in the crystal field with I = 0.2 and with sufficiently 
low temperature after pulse is over. The evolution of vacancies show a very slow (virtually none) 
dynamics well after the electric pulse has been applied since the crystal field barriers are high with 
respect to the thermal noise. 

between different phases is very slow (logarithmic-time dynamics). We observe similar be- 

lavior in our simulations. Indeed, as observed in our earlier MDs without the crystal field 
1 

2|, the vacancy distribution front keeps moving towards the top of the sample gradually 
forming filamentary clusters even after switching the voltage pulse off. If the vacancies 
continue to move for long times after the end of the pulse, it would be difficult to build a 
non- volatile device. These earlier simulations might be relevant for other systems, like Li + , 
Na + or even Ag + ions in Silicon or Germanium. In these systems, it is known that the 
ions continue to move after the pulse is over. Remarkably, we observe in the present MDs 
that the crystal field stabilizes the clusters after they are formed by lowering the mobility 
of vacancies. These non-equilibrium metastable vacancy states relax extremely slowly, since 
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they have to overcome crystal field barriers to relax to states having lower energy. Thus, 
vacancies in memristors can be frozen into a metastable phase for an extremely long time, 
in analogy to vortices staying in non-equilibrium critical states for many years. As a result, 
we see a non- volatile behavior where the vacancies freeze right after the pulse if the crystal 
field is large enough, as illustrated in Fig. (j3J) representing the vacancy distributions on a 
long time scale. 

We conclude that manipulating the vacancy dynamics with the crystal field and the 
vacancy- vacancy interactions could provide a rather rich distribution of phases. We observe 
a tendency towards crystallization of the vacancies, in agreement with the experimental 
devices, where the Ti02 material in the switching channels is crystalline jsj. Our MD 
simulations indicate that competing with the randomizing stochastic forces, some crystal 
field potentials provide the orientational ordering and stabilization of electrically conducting 
channels in memristors that are vital for their nonvolatile memory property. Further 3D MDs 
of real memristors with the number of oxygen vacancies as large as 1000 and above should 
provide a powerful approach for understanding and improving these devices. 
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